function [stt_o,eta_o,smooth_o,ptt_o,state_dr,eta_dr,yfor_o,yferr_o]=kalman_dr1(y,xdata,Z,Cobs,T,R,Qmat,A,addin)
% ==================================================================
% [stt_o,eta_o,smooth_o,ptt_o,state_dr,eta_dr,yfor_o,yferr_o]=kalman_dr1(y,
% xdata,Z,Cobs,T,R,Qmat,A,addin)
% KALMAN_DR1
% 
% This function takes as inputs the matrices of a state space model with no
% error in the observation equation and the Kalman Gain, and F^(-1) matrices from 
% a forward run of the Kalman Filter in order to simulate a single set of draws
% from the state disturbances and the states using a disturbance smoother. 
% It also returns the filetered and smoothed states. 
% It works both with ESTDGSE and SSM type models. 
% Only for NON_DIFFUSE state space with constant variances 
%
% The models is 
% y(t) = Z*s(t) 
% s(t) = T*s(t) + A*xdata(t) + R*n(t)  V( n(t) ) = Q=Qmat'*Qmat   
%
% Where y is [NZ  1], xdata(t) [NE  1]; s(t)[NS 1] and n(t) [NX x 1] 
% 
% Related files: 
% KALMAN_DRALL.m  Generate multiple draws 
% Inputs 
% ======
% y     nobs x NZ stacked vector of observations, nobs is the total number of observations  
% 
% All matrices but Q(t) are assumed to be time invariant to simplify the computation 
% The correct size of input Qmat then depends on whether Q is time varyong or  not. 
% not time varying    Qmat is NX x NX 
%      time varying   Qmat is NX x NX x T where
% 
% Alejandro Justiniano and Giorgio Primiceri (c) April 2005 
% =============================================================

% 1. Initial check for dimension
% -------------------------------
[nobs,nz]=size(y);
[ny,nx] = size(R);
ns=size(T,1);
Q=Qmat'*Qmat;
RQ=R*(Qmat'); 

% 2. Removal of Auxiliary data and means
% ydem
% ---------------------------------------
if ~isempty(xdata);
    %flag_xdata=1;
    if size(xdata,1)~=nobs;error('X and Y must have same number of rows');end
    if size(xdata,2)~=size(A,2);error('Cols A must match Cols of X');end;
    if size(A,1)~=size(T,1);error('Rows A must match dimension state space');end;
    ydem=y-xdata*(A');
else
    %flag_xdata=0;
    ydem=y;
end

if any(Cobs~=0)==1;
    ydem=ydem-repmat((Z*Cobs)',[nobs 1]);
end;

% 3. Initialization
% ==================
% I. Non-diffuse
% PSHAT P0|0
% SHAT  s0|0 
if addin.flag_pzero==1
    pzero=addin.pzero;
else
    pzero=lyapunov_symm(T,RQ*(RQ'));
    pzero=0.5*(pzero+pzero'); 
    if ~isempty(find(isnan(pzero))~=0)
        error('PSHAT contains Nan entries')
    end
end
if addin.flag_szero==1;
    szero=addin.szero;
    if length(szero)~=ns;error('SHAT must be NSx1');end
else
    szero=zeros(ns,1);
end

% 4. Filter YDEM 
% Resulting matrices have _0 for orignal 
% =====================================================

[stt_o,ptt_o,vstar_o,finmat_o,kpartg_o,yfor_o]=kf_forward_sub(ydem',szero,pzero); 
   
% 5. Smooth YDEM 
% Resulting matrices have _0 for orignal 
% ========================================
Ztr=Z'; %clear Z;
Rtr=R'; %clear R;
[eta_o,smooth_o]=kf_backward_sub(ydem',szero,pzero,finmat_o,vstar_o,kpartg_o,stt_o,ptt_o); 

stt_o = stt_o'; 
yferr_o=ydem - yfor_o; 

% ============================================= 
% 6.    Simulate the innovations and the data
% w+
% y+
% s+    alpha+
% =============================================
wplus=randn(nobs,nx)*Qmat;
% This is y+ and alpha+ in DK.
% SZERO_PLUS: first state drawn from a Normal (shat,pshat)
[yplus,stplus,szero_plus]=simssmodel(wplus,T,R,Z,szero,pzero);


% =============================================
% 7.  Forward and backward with YPLUS 
% =============================================
[stt_plus,ptt_plus,vstar_plus,finmat_plus,kpartg_plus]=kf_forward_sub(yplus',szero,pzero); 
[eta_plus,smooth_plus]=kf_backward_sub(yplus',szero,pzero,finmat_plus,vstar_plus,kpartg_plus,stt_plus,ptt_plus); 

% =============================================
% 8. DRAW 
% ==============================================
eta_dr=eta_o   -eta_plus   +wplus;
state_dr=smooth_o-smooth_plus+stplus;

% =========================================================================
% =========================================================================
% =========================================================================

    % =====================================================================
    % Soubroutine 1 
    % KF_FORWARD_SUB 
    % Same as KF_DK but in subroutine mode 
    % STT and PTT DO NOT have the inital S0 and P0 in the first entry 
    % =====================================================================
    function [stt,ptt,vstar,finmat,kpartg,yfor]=kf_forward_sub(yy,a_zero,p_zero)
        vstar=zeros(nz,nobs);
        finmat=zeros(nz,nz,nobs);
        kpartg=zeros(ns,nz,nobs);
        logLnc=zeros(nobs,1);
        yfor=zeros(nz,nobs);
        % Matrices with one additional entry (initialization)
        stt=zeros(ns,nobs+1);
        ptt=zeros(ns,ns,nobs+1);
        stt(:,1)  =a_zero;
        ptt(:,:,1)=p_zero;
        for ii=1:nobs;
            yfor(:,ii)=Z*T*stt(:,ii);
            [stt(:,ii+1),ptt(:,:,ii+1),logLnc(ii),vstar(:,ii),finmat(:,:,ii),kpartg(:,:,ii)]...
                =kf_dk(yy(:,ii),Z,stt(:,ii),ptt(:,:,ii),T,RQ);
        end
        stt=stt(:,2:end); 
        ptt=ptt(:,:,2:end); 
        yfor=yfor';        
    end

% =================================================================
% Subroutine 2
% KF_BACKWARD_SUB
% Backward routine for STATE and disturbance smoother
% Using the inputs of KF_FORWARD_SUB
% =================================================================

    function [etamat,smooth_st]=kf_backward_sub(yy,a_zero,p_zero,finmat,vstar,kpartg,stt,ptt)
          
        etamat   =zeros(nx,nobs);
        smooth_st=zeros(ns,nobs);
        
        rstar=zeros(ns,1);
        rmat=zeros(ny,nobs);
        
        [rstar,etamat(:,end)]=smoothdis(rstar,Q,Rtr,Ztr,finmat(:,:,end),zeros(ns),vstar(:,end));
        smooth_st(:,nobs)=stt(:,end)+ptt(:,:,end)*rstar;
        rmat(:,nobs)=rstar;
        
        for ii=(nobs-1):-1:1;
            [rstar,etamat(:,ii)]=smoothdis(rstar,Q,Rtr,...
                Ztr,finmat(:,:,ii),((T-T*kpartg(:,:,ii)*Z)'),vstar(:,ii));
            smooth_st(:,ii)=stt(:,ii)+ptt(:,:,ii)*rstar;
            rmat(:,ii)=rstar;
        end
        
        % =============================================================
        % As a check perform another disturbance smoother and check can
        % recover the observables 
        % ==============================================================        
        smooth_st2=zeros(ns,nobs); 
        smooth_st2(:,1)=T*a_zero+p_zero*rstar; 
        ycheck=zeros(size(yy)); 
        ycheck(:,1)=Z*smooth_st2(:,1); 
        %R=Rtr'; 
        for ii=2:nobs;
            smooth_st2(:,ii)=T*smooth_st2(:,ii-1) + R*(etamat(:,ii));
            ycheck(:,ii)=Z*smooth_st2(:,ii);
        end
        ydiscrep=max((abs(ycheck'-yy'))); 
        if max(ydiscrep > 1e-4); 
            error('SMOOTH innovations do not recover observables'); 
        end 
        
        % =============================================
        % Check also difference between smooth states 
        if max(max(abs( smooth_st2-smooth_st ))) > 1e-4 
            error('SMOOTH states do not match') 
        end 
   
        etamat=etamat';
        smooth_st=smooth_st';
                    
    end

end 